Lab 7
Shone George Kutty Renjan
2024-10-14
Tasks
Task 1: Set working directory
getwd()## [1] "/Users/shonerenjan/Desktop/ASM/Labs/Lab 7"Task 2: Chi-Square-Statistic
mychisim<-function(n1=10,sigma1=3,mean1=5,iter=1000,ymax=0.1,x=20, y=0.1)
{ # adjust ymax to make graph fit
y1=rnorm(n1*iter,mean=mean1,sd=sigma1)# generate iter samples of size n1
data1.mat=matrix(y1,nrow=n1,ncol=iter,byrow=TRUE) # Each column is a sample size n1
ssq1=apply(data1.mat,2,var) # ssq1 is s squared
w=(n1-1)*ssq1/sigma1^2 #chi-sq stat
hist(w,freq=FALSE, ylim=c(0,ymax), # Histogram with annotation
main=substitute(paste("Sample size = ",n[1]," = ",n1," statistic = ",chi^2)),
xlab=expression(paste(chi^2, "Statistic",sep=" ")), las=1)
lines(density(w),col="Blue",lwd=3) # add a density plot
curve(dchisq(x,n1-1),add=TRUE,col="Red",lty=2,lwd=3) # add a theoretical curve
title=expression(chi^2==frac((n[1]-1)*s^2,sigma^2)) #mathematical annotation -see ?plotmath
legend(x,y,c("Simulated","Theoretical"),col=c("Blue","Red"),lwd=4,lty=1:2,bty="n",title=title) # Legend #
return(list(w=w,summary=summary(w),sd=sd(w),fun="Chi-sq")) # some output to use if needed
}\[ n_1=10, iter=1000,μ_1=10,σ_1=4\]
chisim1=mychisim(n1=10, sigma1=4, mean1=10, iter=1000, y=0.06)
\[n_1=20,iter=1000,μ_1=10,σ_1=4\]
chisim2 = mychisim(n1=20, sigma1 = 4, mean1 = 10, iter = 1000, x=30, y=0.05 )
\[n_1=100,iter=1000,μ_1=10,σ_1=4\]
chisim3=mychisim(n1=100,sigma1=4,mean1=10,iter=1000,ymax=0.06, x=120,y=0.05)
\[ n_1=200,iter=1000,μ_1=10,σ_1=4\]
chisim4=mychisim(n1=200, sigma1=4,mean1=10,iter=1000,ymax=0.05,x=230,y=0.04)
\[n_1=10,iter=1500,μ_1=20,σ_1=10 \]
chisq=mychisim(n1=10, sigma1=10,mean1=20,iter=1500)
hist(chisq$w, col=rainbow(12 ), xlab=expression(paste(chi^2,"Value",sep=" "),las=1))
Task 3: T-Statistic
myTsim<-function(n1=10,sigma1=3,mean1=5,iter=1000,ymax=0.5,x=2,y=0.25){ # adjust ymax to make graph fit
y1=rnorm(n1*iter,mean=mean1,sd=sigma1)# generate iter samples of size n1
data1.mat=matrix(y1,nrow=n1,ncol=iter,byrow=TRUE) # Each column is a sample size n1
sd1=apply(data1.mat,2,sd) # sd
ybar=apply(data1.mat,2,mean) # mean
w=(ybar-mean1)/(sd1/sqrt(n1)) #T stat
hist(w,freq=FALSE, ylim=c(0,ymax), # Histogram with annotation
main=substitute(paste("Sample size = ",n[1]," = ",n1," statistic = ",T," iterations= ",iter)),
xlab=expression(paste(T, "Statistic",sep=" ")), las=1)
lines(density(w),col="Blue",lwd=3) # add a density plot
curve(dt(x,n1-1),add=TRUE,col="Red",lty=2,lwd=3) # add a theoretical curve
title=expression(T==frac((bar(y)-mu),s/sqrt(n1))) #mathematical annotation -see ?plotmath
legend(x,y,c("Simulated","Theoretical"),col=c("Blue","Red"),lwd=4,lty=1:2,bty="n",title=title) # Legend #
return(list(w=w,summary=summary(w),sd=sd(w),fun="T")) # some output to use if needed
}\[ n_1=10,iter=1000,μ_1=10,σ_1=4\]
TS1=myTsim(n=10, sigma1=4, mean1=10, iter=1000)
\[ n_1=20,iter=1000,μ_1=10,σ_1=4\]
TS2=myTsim(n1=20, sigma1=4, mean1=10, iter=1000,ymax=0.50, x=-3.75, y=0.4)
\[ n_1=100,iter=1000,μ_1=10,σ_1=4\]
TS3=myTsim(n1=100, sigma1=4, mean1=10, iter=1000,ymax=0.50, x=-3.75, y=0.4)
\[ n_1=200,iter=1000,μ_1=10,σ_1=4\]
TS4=myTsim(n1=200, sigma1=4, mean1=10, iter=1000,ymax=0.50, x=-3, y=0.4)
\[n_1=10,iter=1500,μ_1=20,σ_1=10\]
T=myTsim(n1=10,mean1 = 20,sigma1 = 10,iter = 1500,x=-5,y=0.3)
hist(T$w, col=rainbow(12), xlab = expression(paste(T, " Value",sep+" "), las=1))
Task 4 Chi-Square Statistic 2
mychisim2<-function(n1=10,n2=14,sigma1=3,sigma2=3,mean1=5,mean2=10,iter=1000,ymax=0.07,x=2,y=0.05,...){ # adjust ymax to make graph fit
y1=rnorm(n1*iter,mean=mean1,sd=sigma1)# generate iter samples of size n1
y2=rnorm(n2*iter,mean=mean2,sd=sigma2)
data1.mat=matrix(y1,nrow=n1,ncol=iter,byrow=TRUE) # Each column is a sample size n1
data2.mat=matrix(y2,nrow=n2,ncol=iter,byrow=TRUE)
ssq1=apply(data1.mat,2,var) # ssq1 is s squared
ssq2=apply(data2.mat,2,var)
spsq=((n1-1)*ssq1 + (n2-1)*ssq2)/(n1+n2-2) # pooled s squared
w=(n1+n2-2)*spsq/(sigma1^2)#sigma1=sigma2, Chi square stat
hist(w,freq=FALSE, ylim=c(0,ymax), # Histogram with annotation
main=substitute(paste("Sample size = ",n[1]+n[2]," = ",n1+n2," statistic = ",chi^2)),
xlab=expression(paste(chi^2, "Statistic",sep=" ")), las=1)
lines(density(w),col="Blue",lwd=3) # add a density plot
curve(dchisq(x,n1+n2-2),add=TRUE,col="Red",lty=2,lwd=3) # add a theoretical curve
title=expression(chi^2==frac((n[1]+n[2]-2)*S[p]^2,sigma^2)) #mathematical annotation -see ?plotmath
legend(x,y,c("Simulated","Theoretical"),col=c("Blue","Red"),lwd=4,lty=1:2,bty="n",title=title) # Legend #
return(list(w=w,summary=summary(w),sd=sd(w),fun="Chi-sq")) # some output to use if needed
}\[ n_1=10,n_2=10,μ_1=5,μ_2=10,σ_1=σ_2=4,iter=1000\]
chisq1=mychisim2(n1=10,n2=10, mean1 = 5,mean2 = 10,sigma1 = 4,sigma2 = 4,iter=1000,x=30)
\[ n_1=20,n_2=10,μ_1=3,μ_2=5,σ_1=σ_2=10,iter=1000\]
chisq2=mychisim2(n1=20,n2=10, mean1 = 3,mean2 = 5,sigma1 = 10,sigma2 = 10,iter=1000,x=40)
\[ n_1=50,n_2=50,μ_1=5,μ_2=10,σ_1=σ_2=4,iter=10000\]
chisq3=mychisim2(n1=50,n2=50, mean1 = 5,mean2 = 10,sigma1 = 4,sigma2 = 4,iter=10000,ymax=0.05,x=130,y=0.04)
\[ n_1=80,n_2=50,μ_1=3,μ_2=5,σ_1=σ_2=10,iter=10000\]
chisq4=mychisim2(n1=80,n2=50, mean1 = 3,mean2 = 5,sigma1 = 10,sigma2 = 10,iter=10000,ymax=0.05,x=160,y=0.04)
\[Default\]
defautchisq= mychisim(iter=10000,x=35,y=0.05)
hist(defautchisq$w,col=rainbow(12), xlab=expression(paste(chi^2,"Value",sep = " "),las=1))
Task 5:
student’s T statistic: \[T = \frac{(\bar{Y}_1 - \bar{Y}_2) - (\mu_1 - \mu_2)}{\sqrt{S_p^2 (\frac{1}{n_1} + \frac{1}{n_2})}}\] \[S^2_p = \frac{(n_1 - 1)S^2_1 + (n_2 - 1)S^2_2}{n_1 + n_2 - 2}\] •̄ Yᵢ represents the sample mean of group i • μᵢ stands for the population mean of group i • nᵢ denotes the sample size from group i • s²ᵢ indicates the sample variance of group i • s²ₚ refers to the pooled variance across all sampless the samples.
myTsim2<-function(n1=10,n2=14,sigma1=3,sigma2=3,mean1=5,mean2=10,iter=1000,ymax=0.5,x=2,y=0.4,...){
y1=rnorm(n1*iter,mean=mean1,sd=sigma1)# generate iter samples of size n1
y2=rnorm(n2*iter,mean=mean2,sd=sigma2)
data1.mat=matrix(y1,nrow=n1,ncol=iter,byrow=TRUE) # Each column is a sample size n1
data2.mat=matrix(y2,nrow=n2,ncol=iter,byrow=TRUE)
ssq1=apply(data1.mat,2,var) # ssq1 is s squared
ybar1= apply(data1.mat,2,mean)
ssq2=apply(data2.mat,2,var)
ybar2=apply(data2.mat,2,mean)
spsq=((n1-1)*ssq1 + (n2-1)*ssq2)/(n1+n2-2) # pooled s squared
w=((ybar1-ybar2)-(mean1-mean2))/sqrt(spsq*(1/n1+1/n2))#sigma1=sigma2, Chi square stat
hist(w,freq=FALSE, ylim=c(0,ymax), # Histogram with annotation
main=substitute(paste("Sample size = ",n[1]+n[2]," = ",n1+n2," statistic = ",T)),
xlab=paste(" T Statistic",sep=""), las=1)
lines(density(w),col="Blue",lwd=3) # add a density plot
curve(dt(x,n1+n2-2),add=TRUE,col="Red",lty=2,lwd=3) # add a theoretical curve
title=expression(T==frac((bar(Y)[1]-bar(Y)[2])-(mu[1]-mu[2]),S[p]*sqrt(frac(1,n[1])+frac(1,n[2])))) #mathematical annotation -see ?plotmath
legend(x,y,c("Simulated","Theoretical"),col=c("Blue","Red"),lwd=4,lty=1:2,bty="n",title=title)# Legend #
return(list(w=w,summary=summary(w),sdw=sd(w),fun="T")) # some output to use if needed
}\[ n_1=10,n_2=10,μ_1=5,μ_2=10,σ_1=σ_2=4,iter=1000\]
doubleT1=myTsim2(n1=10, n2=10,sigma1=4,sigma2=4,mean1=5,mean2=10,iter=1000)
\[ n_1=20,n_2=10,μ_1=3,μ_2=5,σ_1=σ_2=10,iter=1000\]
doubleT2=myTsim2(n1=20,n2=10,mean1=3,mean2=5,sigma1=10,sigma2=10,iter = 1000, ymax=0.7,x=-3.5,y=0.5)
\[ n_1=50,n_2=50,μ_1=5,μ_2=10,σ_1=σ_2=4,iter=10000\]
doubleT3=myTsim2(n1=50,n2=50,mean1=5,mean2=10,sigma1=4,sigma2=4,iter = 10000, ymax=0.7,x=-4.25,y=0.5)
\[ n_1=80,n_2=50,μ_1=3,μ_2=5,σ_1=σ_2=10,iter=10000\]
doubleT4=myTsim2(n1=80,n2=50,mean1=3,mean2=5,sigma1=10,sigma2=10,iter = 10000, ymax=0.7,x=-4,y=0.5)
\[Default=iter=10000\]
Doubletd=myTsim2(iter = 10000,x=-4.0)
hist(Doubletd$w,col = rainbow(10),xlab = expression(paste(T," Value", sep=" "),las=1))
Task 6: The F Statistic
The statistic that the function will calculate: \[𝐹=(S^2_1/𝑆^2_2)(𝜎^2_2/𝜎^2_1)\] No assumptions
myFsim2<-function(n1=10,n2=14,sigma1=3,sigma2=2,mean1=5,mean2=10,iter=1000,ymax=0.9,x=6,y=0.5,...){
y1=rnorm(n1*iter,mean=mean1,sd=sigma1)# generate iter samples of size n1
y2=rnorm(n2*iter,mean=mean2,sd=sigma2)
data1.mat=matrix(y1,nrow=n1,ncol=iter,byrow=TRUE) # Each column is a sample size n1
data2.mat=matrix(y2,nrow=n2,ncol=iter,byrow=TRUE)
ssq1=apply(data1.mat,2,var) # ssq1 is s squared
ssq2=apply(data2.mat,2,var)
#spsq=((n1-1)*ssq1 + (n2-1)*ssq2)/(n1+n2-2) # pooled s squared
w=ssq1*sigma2^2/(ssq2*sigma1^2) #
hist(w,freq=FALSE, ylim=c(0,ymax), # Histogram with annotation
main=substitute(paste("Sample size = ",n[1]+n[2]," = ",n1+n2," statistic = ",F)),
xlab=paste("F Statistic",sep=""), las=1)
lines(density(w),col="Blue",lwd=3) # add a density plot
curve(df(x,n1-1,n2-1),xlim=c(0,6),add=TRUE,col="Red",lty=2,lwd=3) # add a theoretical curve
title=expression(F==frac(s[1]^2,s[2]^2)*frac(sigma[2]^2,sigma[1]^2)) #mathematical annotation -see ?plotmath
legend(x,y,c("Simulated","Theoretical"),col=c("Blue","Red"),lwd=4,lty=1:2,bty="n",title=title)# Legend #
return(list(w=w,summary=summary(w),sd=sd(w),fun="F")) # some output to use if needed
}$$$$
F1=myFsim2(n1=10, n2=10,sigma1=4,sigma2=4,mean1=5,mean2=10,iter=1000)
\[ n_1=20,n_2=10,μ_1=3,μ_2=5,σ_1=σ_2=10,iter=1000\]
F2=myFsim2(n1=20, n2=10,sigma1=10,sigma2=10,mean1=3,mean2=5,iter=1000)
\[ n_1=50,n_2=50,μ_1=5,μ_2=10,σ_1=σ_2=4,iter=10000\]
F3=myFsim2(n1=50, n2=50,sigma1=4,sigma2=4,mean1=5,mean2=10,iter=10000,ymax=1.5)
\[n_1=80,n_2=50,μ_1=3,μ_2=5,σ_1=σ_2=10,iter=10000\]
F4=myFsim2(n1=80, n2=50,sigma1=10,sigma2=10,mean1=3,mean2=5,iter=10000,ymax=1.7)
Fdv = myFsim2(iter=10000)
hist(Fdv$w, col=rainbow(8))
Task 7: Get File
library(MATH4753SHON2024)
use_data(fire)
knitr::kable(head(use_data(fire)))| DISTANCE | DAMAGE |
|---|---|
| 3.4 | 26.2 |
| 1.8 | 17.8 |
| 4.6 | 31.3 |
| 2.3 | 23.1 |
| 3.1 | 27.5 |
| 5.5 | 36.0 |